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Resumen 

Las partes estructurales de los reactores nucleares constituyen sistemas mecánicos complejos, que pueden 
vibrar en un conjunto de frecuencias propias, cuando sufren una perturbación apropiada. Se producen así 
variaciones cíclicas en el estado de deformación de los materiales, incluyendo perturbaciones en su 
densidad, las cuales pueden afectar la reactividad. Una alteración de la reactividad afecta la potencia 
térmica del reactor, modificando el campo de temperaturas al que se encuentran sometidos los materiales 
mencionados. Si la variación en el campo de temperaturas es lo suficientemente rápida, el acoplamiento 
termo-mecánico puede provocar variaciones rápidas en el estado de deformaciones, el cual, a su vez, 
afecta la reactividad, y así sucesivamente. Este acoplamiento entre las vibraciones mecánicas de la 
estructura y los materiales, con las oscilaciones en la potencia del reactor, no solamente no parece ser 
pasible de ser excluido a priori, sino que parece haber estado presente en alguna etapa durante ciertos 
incidentes o accidentes que acompañaron el desarrollo de la tecnología de reactores. El propósito de esta 
comunicación es doble. (a) Revisar y generalizar algunos modelos matemáticos propuestos para describir 
el acoplamiento núcleo-mecánico. (b) Discutir algunas condiciones en las que podrían surgir 
inestabilidades significativas, incluyendo oscilaciones de gran amplitud en la potencia, acompañadas por 
vibraciones mecánicas lo bastante pequeñas como para no resultar excluidas por el diseño mecánico 
convencional. En vísperas de un resurgimiento de la generación núcleo eléctrica, todo esto podría resultar 
de interés para discutir ciertos aspectos de la seguridad física de los reactores en los que los opositores a 
la tecnología de reactores nucleares están poniendo énfasis. 


(A) Introducción: 


El presente trabajo tiene por objetivo discutir un tipo de inestabilidad oscilatoria en 
la potencia, que puede aparecer en reactores térmicos si se dan ciertas 
circunstancias. Esas circunstancias son muy poco probables, en particular visto el 
estado actual del arte en lo referente al diseño y operación de reactores, pero en 
principio son posibles. 

Al parecer este tipo de inestabilidad oscilatoria habría formado parte de los procesos 
durante ciertos incidentes o accidentes ocurridos tanto en reactores de investigación 
como en reactores de potencia, desde la serie de experimentos con reactores de 
investigación conocidos como BORAX y SPERT realizados en la década del 
cincuenta del siglo pasado, hasta el accidente en el reactor de potencia del tipo 
RBMK de Chernobyl, en 1986. 

No obstante, la mencionada inestabilidad, caracterizada por un acoplamiento entre 
los procesos nucleares y las vibraciones mecánicas del núcleo del reactor mediada a 
través de procesos termo-elásticos, ha sido escasamente estudiada en la literatura 
abierta sobre ciencia y tecnología nucleares. La primera mención pública parece 
haber sido el trabajo de A.S.Thompson, “Study of reactor kinetics”, publicado en las 
actas de la reunión anual de 1962 de la American Society of Mechanical Engineers 
(Thompson y Thompson, 1988). 


La idea básica es bastante simple. Supongamos un reactor en estado estacionario, 
bien diseñado. Si por alguna razón se produce una perturbación que se acompaña de 
un aumento de potencia, entran en juego diversos mecanismos pasivos de 
retroalimentación negativa, inherentes a la física del reactor, unos prácticamente 
instantáneos (como el Doppler) y otros más o menos retardados (como algunos de 
los asociados con el moderador y el fluido refrigerante), que tienden a limitar ese 
incremento sin necesidad de recurrir a los mecanismos activos de control. 

Un incremento en la potencia aumenta la temperatura y produce dilatación y 
distorsión en los materiales del reactor. Cuando la variación en el campo de 
temperatura es lo suficientemente rápida, los efectos inerciales producen un retardo 
entre el incremento de temperatura y la respuesta mecánica, con lo cual el efecto 
sobre la reactividad debida a la respuesta mecánica en los materiales también se 
produce con un retardo. Aun cuando la distorsión de los materiales del reactor en 
respuesta a un incremento de temperatura, por sí misma en condiciones estáticas 
tienda a disminuir la reactividad, debido al retardo cabe la posibilidad de que en 
ciertas condiciones se exciten oscilaciones en la potencia, acopladas con las 
vibraciones mecánicas de los materiales y asociadas a los mencionados efectos 
inerciales. 


Si los parámetros que caracterizan los subsistemas mecánico y neutrónico se 
encuentran dentro de un intervalo de valores adecuado, en particular si la disipación 
de energía mecánica durante las vibraciones no es suficiente y la potencia de estado 
estacionario es lo bastante grande, cabría esperar que como consecuencia de una 
perturbación la potencia pueda presentar oscilaciones muy rápidas e inestables (con 
amplitud creciente), como las que Stephenson (1965) menciona que fueron 
observadas durante algunos de los experimentos SPERT y como las que predicen 
algunos de los modelos matemáticos que contemplan explícitamente este 
acoplamiento entre fenómenos termomecánicos y nucleares. 


Los modelos analíticos clásicos, a parámetros concentrados, que se desarrollaron 
para describir accidentes de reactividad, como el de Nordheim-Fuchs, el de Nyer y 
sus diferentes variantes (Weinberg y Wigner, 1962; Smets, 1962; Hetrick, 1971), 
operan en la misma escala de tiempo (décimas de segundo) en la que se produce el 
acoplamiento núcleo-mecánico pero no permiten abarcar los aspectos mecánicos del 
fenómeno. Solo el modelo analítico de Bethe y Tait, desarrollado para describir un 
accidente de reactividad muy severo en un reactor de neutrones rápidos (Ram, 
1983), abarca algunos de los aspectos mecánicos necesarios. Por ello y por haber 
sido construido a parámetros distribuidos, constituye un antecedente valioso para 
utilizar en un enfoque mediante métodos de análisis modal no lineal, pese a que su 
formulación es tal que no permite describir oscilaciones mecánicas. 

Los modelos utilizados para describir las inestabilidades termo-hidraúlicas en los 
reactores térmicos de potencia, incluyendo las oscilaciones debidas al “flashing”, a 
las alternancias entre patrones de flujo, las de Ledinneg, las dinámicas y las térmicas 
(Tyror y Vaughan, 1970; Rust, 1983) operan por lo general en escalas de tiempo 
mucho mayores que la que aquí interesa. 


El propósito de este trabajo es doble: 

(1) Revisar dos modelos analíticos clásicos, a parámetros concentrados (cinética 
puntual) que se pueden utilizar para estudiar inestabilidades en la potencia del 
reactor. 


(2) Generalizar dichos modelos para poder discutir algunos aspectos del 
acoplamiento núcleo-mecánico y sugerir la aplicación de modelos basados en 
métodos de análisis modal no lineal para modelos a parámetros distribuidos del 
acoplamiento núcleo-mecánico. 


Puesto que los reactores de potencia requieren un abordaje más detallado, 
considerando cada tipo de reactor por separado y recurriendo a códigos de cálculo 
numérico que no siempre pueden abarcar el tipo de fenómeno que nos ocupa (van 
Dam, 1992), un estudio preliminar como el sugerido podría resultar de utilidad para 
diseñar estudios más profundos y realistas, utilizando códigos de cálculo adecuados. 


No obstante, deben tenerse presentes las limitaciones de los modelos de cinética 
puntual para representar procesos en escalas de tiempo muy cortas, tanto en 
reactores de investigación o producción como en reactores de potencia. 


(B) Un modelo de cinética puntual con retardo 


Podemos comenzar modificando un modelo muy simple que Akcasu, Lellouche y 
Shotkin (1971, p.441) utilizaron para analizar la posibilidad de tiempos de escape 
finitos en la potencia del reactor. 

Supongamos que hasta el instante inicial el reactor opera a la potencia estacionaria 
P, . En ese instante sufre una perturbación instantánea que lleva la potencia a P(0). 
Como nos interesa estudiar la posibilidad de inestabilidades oscilatorias en escalas 
de décimas de segundo, podemos suponer que la potencia viene dada por la 
ecuación: 
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Aquí el término Er, representa la potencia (constante en esa escala) asociada a 


los neutrones retardados. Como es usual, Pf representa la fracción de neutrones 


retardados y A representa el tiempo promedio entre generaciones de neutrones. 
Suponemos ahora que la reactividad presenta un mecanismo de retroalimentación 
instantáneo y un mecanismo retardado, siendo 7 el retardo: 


p=6p, +0, (PB) +0, (PU=T).UU=r) [2] 


La modificación de reactividad óp, debida a causas externas, como movimientos 
en las barras de control u otros eventos, se considerará nula o se supondrá que se 
anula luego de producir una perturbación en la potencia del reactor que la aleja de su 
valor de régimen. La función U(t—7) es la función escalón unidad de Heavside. 
Indica que el mecanismo retardado comienza a actuar luego de transcurrido un 
tiempo 7. Entonces P(t)=P, es solución para t<0. Los coeficientes de 
retroalimentación instantánea y retardada de potencia se representan por los 
parámetros a, y «a, respectivamente. En todos los casos suponemos que «, es 
negativo (lo cual permite tener en cuenta la retroalimentación negativa originada en 
el efecto Doppler en el combustible). 
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aproximación lineal en torno a ¿=0 , obtenemos a partir de [1] y [2]: 


Introduciendo la nueva variable de estado ¿(t)= y efectuando una 


E a dla, P,- B)éli)+a,.P, ¿(ir JU (1) 13] 


poa : dxki 
Esta es una ecuación diferencial con retardo del tipo Ext) = Axlt)+ B.x(t—7), que 
Í 


fue bien estudiado por Hayes (1950). En lo que sigue utilizaremos algunos de sus 

resultados aplicándolos al estudio de la estabilidad de un reactor nuclear. 

Las soluciones de este tipo de ecuaciones con retardo se construyen combinando 

linealmente términos exponenciales e”**, donde A es raíz de la ecuación 4 = 4+B.e”* 

a Py =P 
A 


a,.P, E 
. En el caso del reactorx=¿ , 4= y B= o Como «, es negativo el 


coeficiente del término no retardado es siempre negativo. La estabilidad del reactor 
depende de la raíz cuya parte real es máxima, y las oscilaciones se producen si la parte 
imaginaria es distinta de cero. 

A A 
Si AN ,P(t)regresa a P, , sin oscilar, luego de la perturbación. Pero si AN la 
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La función arco coseno se toma entre 0 y 7. Si el retardo supera ese valor crítico 
T.(que decrece con el incremento de la potencia estacionaria P, ), la potencia P oscila 


[4]. 


potencia regresa a su valor estacionario solamente si T<T,= 


en torno a P,con amplitud creciente hasta que la aproximación lineal deja de ser 


aplicable. En las proximidades de la transición a la inestabilidad es posible estimar el 
período de oscilación 7, de la potencia: 27 ,<T, <4r,si Ay B son ambos negativos, y 


27, <T, si B es positivo. 


Cuando el retardo toma su valor crítico 7, la frecuencia de oscilación de la potencia 
: 1 
viene dada por: O, =VB?*- 4? = alla, PY —La,|.., + PB) [5] 


Aún cuando ambos coeficientes de retroalimentación de la potencia sean negativos, la 
potencia del reactor se inestabiliza oscilando con amplitud creciente cuando el retardo 


el reactor es estable a 


, resulta que si [a >|a, 


supera 7.. Como ! 
ce 


> 


<la, 


cualquier potencia estacionaria de operación, según el modelo lineal. Pero si a, 


existe una potencia umbral P, = P,,, por debajo de la cual el reactor es estable, pero por 
encima se inestabiliza oscilando con amplitud creciente si el retardo es mayor que el 


ÓS 


correspondiente retardo crítico 7, : 


OS 
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(C) Modificación del modelo del modelo de Wigner-Weinberg-Thompson 


El modelo simple que construimos muestra cómo pueden surgir inestabilidades 
oscilatorias pero no plantea mecanismos explicativos. El siguiente paso será, entonces, 
refinar el modelo para exhibir en forma explícita el acoplamiento núcleo-mecánico. 
Siguiendo a Thompson (Thompson y Thompson, 1988) cuando modificó un modelo de 
cinética puntual de Wigner y Weinberg (Weinberg y Wigner, 1962), vamos a suponer 
que la reactividad se puede expresar en términos de una temperatura  T asignada 
globalmente al núcleo del reactor y de una variable y que representa la variación de 
una densidad global de los materiales del núcleo del reactor, relativa a su valor en el 
estado de régimen a potencia P, . 

El efecto sobre la reactividad de una variación en la temperatura se describe mediante el 
coeficiente q, y el efecto debido a la variación en la densidad mediante a ,. 

Ahora generalizaremos el modelo de Wigner-Weinberg-Thompson, añadiendo un 
término dependiente del valor instantáneo de la potencia, para tener en cuenta los 
mecanismos de retroalimentación que operan a nivel del combustible y que se pueden 
considerar instantáneos, como el Doppler. Cabe esperar que esos mecanismos no se 
puedan relacionar en forma directa con la temperatura 7 que se asigna al núcleo del 
reactor considerado globalmente. Introduciendo entonces un coeficiente de 
realimentación instantánea de potencia, a, , la reactividad puede escribirse así: 

p=0p,+0 (P-P,)+ a, (T-T,)+a,.y [7] 


Igual que antes, supondremos que óp, se puede considerar nula. 
4 ; dT 
La potencia evoluciona según la [1] y la temperatura según  C. en =P-P, [8] 
É 


Esta relación expresa que todo exceso de potencia térmica por encima de P, no puede 


ser evacuado por el refrigerante del núcleo y se acumula entonces en los materiales con 
una velocidad de aumento de la temperatura inversamente proporcional a la capacidad 
calorífica C del núcleo. Por el contrario, toda disminución de potencia por debajo de 
P, se acompaña de una disminución de 7. 


La variación de densidad y la temperatura se supone que se conectan por la relación: 


2 
PE ES [9] 


En esta ecuación c,, es un coeficiente de fricción mecánica, «,,es una frecuencia natural 


de vibración y bes un coeficiente que describe globalmente el efecto termo-elástico en 
los materiales del reactor. 

Comenzando con la [1], teniendo en cuenta la [7] y utilizando la [8] para eliminar la 
temperatura se obtiene, luego de varios pasos de cálculo: 


d? le 


z E Pd(P-P P(P-BY a 
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Aquí, por definición v= E es la velocidad de variación de la densidad relativa. 
14 


Eliminando la temperatura de la [9], se obtiene 


[11] 


a EE v= 
dt dt 8 


d?y dv >, 0, DP, (5 ia 
Las ecuaciones [10] y [11] describen el acoplamiento núcleo-mecánico. La [10] es una 
generalización de la ecuación correspondiente del modelo de Thompson, para tener en 
cuenta los mecanismos instantáneos de retroalimentación negativa. 


; : ; : P sd 
Es conveniente introducir la nueva variable x = me) . La ecuación [10] se puede re- 
0 


») 
197 
escribir entonces de esta forma: gor + f(x) EA +0, g(x)==v [12] 
dt dt A 
Por definición:  f(x)=c,.e*+c,.e* y g(x)=e*" 1 
a, |.P, = Ay P, 
Los parámetros C, Cp, Op. se definen así: c,= E q == a = E 
A A AC 


El coeficiente «, se supone es negativo (retroalimentación por efecto de temperatura 
negativa) mientras que el coeficiente «, se toma positivo (retroalimentación por efecto 
de vacío también negativa). La [12] es la ecuación de un oscilador nuclear no lineal 


amortiguado y acoplado con el oscilador mecánico [11], cuya ecuación se puede re- 
escribir así: 
2 
d*y dv ) o, b.P, 
Ke, +0, v= 2 gl) 


13 
de "dt C 115] 


Si se efectúa una aproximación lineal en torno ax= 0, se obtienen las ecuaciones: 


de o 2, 4, dy dv, O, b.P, 
C +07), +0, O EAN 


: 2 Q,b.P, e 
Introduciendo 0,” = Ae y eliminando v entre esas dos ecuaciones lineales resulta: 


dx 
en 


4 3 2 
a +0, +0) lento +cp)rop +0, e o +lc +0), 


Esta ecuación difiere de una análoga que deducen, por otro camino, Thompson y 
Thompson (1988), solamente en el coeficiente c, , que en la de ellos no aparece. 


Todos los coeficientes del correspondiente polinomio de estabilidad (o característico) 
a, +a,z+a,z" +a,27 +2z* son positivos, por lo cual podemos asegurar la estabilidad 
imponiendo la restricción a,”.a, <a(a,.a, —a,) que se desprende del criterio de 
Routh-Hurwitz (Akcasu y otros, 1971, pp. 255 a 261). 


Los parámetros del oscilador núcleo-mecánico dependen de los materiales y de la 
configuración específica de cada núcleo de reactor. No obstante, si se asumen los 
intervalos de valores para esos parámetros que consideró Thompson, la condición de 
estabilidad se puede simplificar, también en el caso del modelo generalizado que 
presentamos en este trabajo. 


% E . r 2 
Dicha condición se puede expresar aproximadamente así: 0, <C,, le, +C ls) [14] 


oz .b.P, 
Reintroduciendo los parámetros originales la [14] queda así: aa pS (6 + a, 12 le, 


a, b b.a 


> la, .c,, existe un umbral de potencia: B,= a [15] 
1-2 la, 
ba, 


Si 


Si la potencia P, supera ese umbral el reactor se hace inestable. Eso significa que al 


menos una raíz del polinomio característico tiene parte real positiva. Pero como todos 
los coeficientes de dicho polinomio son positivos, esa raíz no puede ser real. Por lo 
tanto la potencia debe aumentar oscilando con amplitud creciente hasta que la 


.C,, , ese umbral de potencia que 


0% b 
aproximación lineal deja de ser aplicable. Si e < pa 


marca la transición a la inestabilidad no existe. 
(D) Discusión y conclusiones 


(1) La condición de estabilidad [14] sugiere que aumentar la fricción en el subsistema 
mecánico y aumentar el coeficiente de retroalimentación instantánea en el subsistema 
nuclear estabiliza el reactor e impide la aparición de las excursiones de potencia 
rápidamente oscilantes con amplitudes crecientes que hemos estado considerando. El 


_le A, 


coeficiente Cc, = Ue que no aparece en el modelo de Thompson del acoplamiento 


núcleo-mecánico, es un parámetro clave en este sentido. Al ser tenido en cuenta se 
podrían modificar significativamente las conclusiones a las que algunas personas han 
llegado en relación con este problema. 

No obstante, una aproximación más precisa exige considerar cada tipo de reactor por 
separado junto con las modificaciones que los parámetros cinéticos sufren con el 
quemado y la gestión del combustible en el núcleo. 


(2) Comparando las fórmulas para la potencia umbral en el modelo con retardo 

(ecuación [6]) y en el modelo del acoplamiento núcleo mecánico (ecuación [15]) se 

desprende que el modelo con retardo puede describir el acoplamiento núcleo-mecánico 
b.a 


si se define el parámetro «, de realimentación retardada así: a == [16] 


r 


ad 
Una vez obtenida esta equivalencia, podemos estimar la frecuencia de oscilación del 
reactor y el tiempo de retardo crítico cuando se alcanza el umbral de inestabilidad. 
Sustituyendo el valor de a, que aparece en la [5] por su expresión como función del 
coeficiente termo-elástico, el coeficiente de dilatación térmica, el coeficiente de fricción 
mecánica y la capacidad calorífica del núcleo tal como aparece en [16]: 


[17] 


Sustituyendo en la [4] se tiene la siguiente estimación del retardo crítico: 


arccos 


(a,L?, + 8) 
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b.a 
T.=A = [18] 
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(3) En principio, los parámetros del modelo de acoplamiento núcleo-mecánico se 
pueden estimar a partir de mediciones obtenidas en circunstancias apropiadas. Los 
datos generados a través de los experimentos SPERT y BORAX se deberían poder 
utilizar para una caracterización de los parámetros de los subsistemas nuclear y 
mecánico en los correspondientes reactores de investigación. Cabe esperar que en 
reactores de investigación, la formación de un número suficiente de burbujas en el agua 
de la piscina en la que se encuentran sumergidas las placas de combustible, pueda 
producir un sistema mecánico en el cual una masa grande de agua en combinación con 
un elemento equivalente a un resorte blando, formado por las burbujas, pueda comenzar 
a vibrar con una frecuencia natural lo suficientemente pequeña como para que el 
acoplamiento núcleo-mecánico sea posible (las frecuencias de oscilación del subsistema 
nuclear desacoplado se encuentran comprendidas entre 0 y1000 Hz, según la estimación 
de Thompson). Este mecanismo permitiría una interpretación concreta del significado 
de los parámetros del subsistema mecánico. En el caso de los reactores de potencia, las 
vibraciones de los sólidos estructurales no deberían dejarse afuera de la consideración, 
al estudiar la posibilidad de acoplamiento núcleo-mecánico. (Posiblemente una parte de 
los procesos asociados con el accidente de Chernobyl tengan relación con esto). 


(4) Puesto que la mayor parte de los códigos de cálculo de reactores no contemplan el 
estudio de pequeñas vibraciones mecánicas en sólidos acopladas con el resto de los 
subsistemas del reactor, podría resultar de cierto interés extender el encare analítico del 
acoplamiento núcleo-mecánico formulando modelos matemáticos relativamente simples 
pero a parámetros distribuídos. Lo más sencillo sería un modelo no lineal, basado en el 
reactor en forma de losa (““slab reactor”), con un solo grupo de neutrones, compuesto 
por una matriz sólida homogénea, imponiendo condiciones de borde al campo de 
temperaturas y a la deformación que se correspondan con el modelo a parámetros 
concentrados. Los problemas de inestabilidad se podrían encarar entonces, al menos en 
forma tentativa, aplicando técnicas de análisis modal no lineal (Suárez-Ántola, 2005ayb; 
Haken, 1983). Los resultados deberían ser de utilidad tanto para el diseño de corridas 
de simulación digital empleando códigos de cálculo como para el diseño de 
experimentos. 


(5) En este trabajo se dedujeron pero no se estudiaron las ecuaciones no lineales [12] y 
[13] que describen el acoplamiento de dos osciladores. Uno lineal, correspondiente al 
subsistema mecánico pero acoplado en forma no lineal al otro subsistema, y el otro no 
lineal correspondiente al subsistema nuclear. 


Thompson y Thompson (1988) efectuaron un estudio numérico de los transitorios del 
modelo de Wigner-Weinberg-Thompson, a partir de la discretización de una ecuación 
diferencial de cuarto orden, que no tiene en cuenta la influencia del parámetro c, que 
se introdujo en el presente trabajo. 


2 
Si se considera el oscilador nuclear desacoplado, e EL (E +0, .2(x)=0 resulta 


que x tiende a cero luego de ser perturbada (P(t) tiende a P,). 

La función de amortiguamiento, siempre positiva, es f (x) =C,e *+cp.e”*. Si se 
despreciac,, el amortiguamiento tiende a cero cuando la potencia aumenta, pero si no 
se lo desprecia el amortiguamiento se recupera e incluso crece exponencialmente con 


P, ; , de é 
(E) Esto sugiere que vale la pena estudiar los transitorios mediante las 
0 


ecuaciones con c, >0. 


(6) Finalmente, un comentario sobre las limitaciones del uso de modelos 
homogeneizados de cinética puntual en la descripción del acoplamiento núcleo-termo- 
mecánico en reactores nucleares. 

El acoplamiento entre el campo homogeneizado de temperaturas y el campo 
homogeneizado de densidades de los materiales se produce, de acuerdo al modelo 
matemático considerado en este trabajo, a través de un efecto termo-elástico 
homogeneizado. 

El cambio de temperaturas que inicia el acoplamiento se produce en el combustible 
nuclear. Si la escala de tiempo de esta variación de temperatura es por lo menos un 
orden de magnitud inferior a la escala de tiempo de transporte de calor desde el 
combustible al refrigerante, no cabe esperar que el acoplamiento núcleo-termo- 
mecánico se produzca, con las consecuencias globales sobre la estabilidad del reactor 
que sugiere el análisis matemático del modelo de cinética puntual. 
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Acoplamiento  núcleo-termo-mecánico en reactores 
nucleares: 

(D) Modelos matemáticos basados en la cinética puntual y 
análisis de la dinámica del reactor mediante una 
aproximación lineal. 


Roberto Suárez Ántola 


Resumen 

Las partes estructurales de los reactores nucleares constituyen sistemas mecánicos 
complejos, que pueden vibrar en un conjunto de frecuencias propias, cuando sufren una 
perturbación apropiada. Se producen así variaciones cíclicas en el estado de deformación 
de los materiales, incluyendo perturbaciones en su densidad, las cuales pueden afectar la 
reactividad. Una alteración de la reactividad afecta la potencia térmica del reactor, 
modificando el campo de temperaturas al que se encuentran sometidos los materiales 
mencionados. Si la variación en el campo de temperaturas es lo suficientemente rápida, el 
acoplamiento termo-mecánico puede provocar variaciones rápidas en el estado de 
deformaciones, el cual, a su vez, afecta la reactividad, y así sucesivamente. Este 
acoplamiento entre las vibraciones mecánicas de la estructura y los materiales, con las 
oscilaciones en la potencia del reactor, no solamente no parece ser pasible de ser excluido 
a priori, sino que parece haber estado presente durante ciertos incidentes o accidentes que 
acompañaron el desarrollo de la tecnología de reactores. Podría explicar en buena medida 
las inestabilidades recientemente observadas en los reactores Maple, destinados a la 
producción de radioisótopos. El propósito de esta primera comunicación es revisar 
algunos modelos matemáticos que pueden aplicarse a una descripción genérica del 
acoplamiento núcleo-termo-mecánico y discutir algunas condiciones en las que podrían 
surgir inestabilidades significativas, incluyendo oscilaciones de gran amplitud en la 
potencia, acompañadas por vibraciones mecánicas lo bastante pequeñas como para no 
resultar excluidas por el diseño mecánico convencional. 


(A) Introducción 


En los reactores nucleares destinados a investigación, producción de 
radioisótopos o conversión de la energía de la fisión nuclear en energía eléctrica, 
puede aparecer un tipo de inestabilidad oscilatoria en la potencia denominada 
inestabilidad por acoplamiento núcleo-termo-mecánico. Las circunstancias en las 
que aparecen estas oscilaciones, rápidas y de amplitud creciente, son poco 
frecuentes. No obstante, en sistemas de la complejidad de un reactor nuclear las 
condiciones que conducen a la desestabilización por acoplamiento núcleo-termo- 
mecánico pueden darse, aún teniendo en cuenta el estado actual del arte en lo que 
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respecta al diseño y la operación de reactores [13], [14]. Debido a las 
consecuencias serias que puede acarrear, la posibilidad de aparición de este tipo 
de inestabilidad no debería descartarse sin efectuar un estudio detallado en cada 
caso concreto. 

La idea básica que permite comprender y prever el acoplamiento núcleo-termo- 
mecánico es bastante simple. Supongamos un reactor en estado estacionario, bien 
diseñado. Si por alguna razón se produce una perturbación que se acompaña de 
un aumento de potencia, entran en juego diversos mecanismos pasivos de 
retroalimentación negativa, inherentes a la física del reactor, unos prácticamente 
instantáneos (como el Doppler) y otros más o menos retardados (como algunos de 
los asociados con el moderador y el fluido refrigerante), que tienden a limitar ese 
incremento en la potencia sin necesidad de recurrir a los mecanismos activos de 
control presentes en todo reactor nuclear [5], [10], [20]. 

Un incremento en la potencia aumenta la temperatura y produce dilatación y 
distorsión en los materiales del reactor. Cuando la variación en el campo de 
temperatura es lo suficientemente rápida, los efectos inerciales producen un 
retardo entre el incremento de temperatura y la respuesta mecánica, con lo cual el 
efecto sobre la reactividad debida a la respuesta mecánica en los materiales 
también se produce con un retardo. Aún cuando la distorsión de los materiales del 
reactor en respuesta a un incremento de temperatura, por sí misma en condiciones 
estáticas tienda a disminuir la reactividad, debido al retardo cabe la posibilidad 
de que en ciertas condiciones se exciten oscilaciones en la potencia, acopladas 
con las vibraciones mecánicas de los materiales y asociadas a los mencionados 
efectos inerciales. Si los parámetros que caracterizan los subsistemas mecánico 
y neutrónico se encuentran dentro de un intervalo de valores adecuado, en 
particular si la disipación de energía mecánica durante las vibraciones no es 
suficiente y la potencia de estado estacionario es lo bastante grande, cabría 
esperar que como consecuencia de una perturbación la potencia pueda presentar 
oscilaciones muy rápidas e inestables, con amplitudes crecientes. En suma, la 
inestabilidad se caracteriza por un acoplamiento entre los procesos nucleares y las 
vibraciones mecánicas del núcleo del reactor mediada a través de procesos termo- 
elásticos. 


Existe evidencia sobre la aparición de este tipo de inestabilidad oscilatoria 
durante ciertos incidentes O accidentes ocurridos tanto en reactores de 
investigación como en reactores de potencia, desde la serie de experimentos con 
reactores de investigación conocidos como BORAX y SPERT realizados en la 
década del cincuenta del siglo pasado [11], hasta el accidente en el reactor de 
potencia del tipo RBMK de Chernobyl, en 1986 [21]. Recientemente, este tipo de 
inestabilidad parece haber sido una de las causas principales de la imposibilidad 
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de poner en funcionamiento los grandes reactores canadienses Maple, cada uno 
de los cuales fue diseñado con el propósito de satisfacer la demanda mundial de 
radioisótopos (de reactor) para aplicaciones médicas e industriales [2]. Si la 
inestabilidad por acoplamiento núcleo-termo-mecánico se llega a manifestar en 
un reactor, puede exigir un rediseño del núcleo (como ocurrió en el caso de los 
reactores Maple, por decisión de la Autoridad Reguladora Nuclear de Canadá). 
Pese a ello, la mencionada inestabilidad ha sido escasamente estudiada en la 
literatura abierta sobre ciencia y tecnología nucleares. 


Los modelos matemáticos que se utilizan corrientemente, tanto en el diseño como 
en el control de la operación de reactores, para describir las inestabilidades 
puramente termo-hidráulicas, o las que resultan de la interacción entre fenómenos 
termo-hidráulicos y neutrónicos, operan en escalas de tiempo mucho mayores que 
la que aquí interesa, por lo que no resultan aplicables. Sobre todo esto se dispone 
una vasta bibliografía, la mayor parte de la cual se relaciona con los reactores de 
agua en ebullición (BWR). Las oscilaciones en la potencia que preocupan a los 
diseñadores y a los operadores de las centrales nucleares basadas en BWR no 
tienen nada que ver con el tipo de problema que nos ocupa, aunque pueda existir 
un cierto parecido, desde el punto de vista formal, en las ecuaciones de los 
modelos simplificados que se emplean en uno y otro caso [6]. 


Los modelos analíticos clásicos, a parámetros concentrados, que se desarrollaron 
para describir accidentes de reactividad, como el de Nordheim-Fuchs, el de Nyer 
y sus diferentes variantes operan en la misma escala de tiempo (décimas de 
segundo) en la que se produce el acoplamiento núcleo-termo-mecánico pero no 
permiten abarcar los aspectos mecánicos del fenómeno [5], [9], [21], [22]. Solo 
el modelo analítico de Bethe y Tait, desarrollado para describir un accidente de 
reactividad muy severo en un reactor de neutrones rápidos abarca algunos de los 
aspectos mecánicos necesarios [7], [10]. Por ello y por haber sido construido a 
parámetros distribuidos, constituye un antecedente valioso, pese a que su 
formulación es tal que no permite describir oscilaciones mecánicas. 


La primera mención pública de la posibilidad de lo que denominamos 
acoplamiento núcleo-termo-mecánico y el primer modelo matemático de dicho 
acoplamiento aparecieron en un trabajo de A.S.Thompson presentado en la 
reunión anual de 1962 de la American Society of Mechanical Engineers [17]'. 


! Bien recibido por los ingenieros mecánicos, este trabajo (de un doctor ingeniero 
mecánico entrenado como ingeniero nuclear a partir de 1946, partícipe activo en el 
diseño y la construcción de las primeras centrales para la conversión núcleo-eléctrica en 
USA, y coautor, con Oliver Rodgers, de uno de los primeros y más estudiados textos de 
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Esta primera comunicación se propone: 

(1) Plantear un modelo analítico inespecífico, basado en una ecuación 
diferencial con retardo, que se puede utilizar para estudiar inestabilidades en la 
potencia del reactor en la escala de tiempo del acoplamiento termo-núcleo- 
mecánico. 

(2) Presentar un modelo reciente que generaliza el modelo original de 
Thompson y que permite deducir expresiones analíticas para los parámetros del 
modelo basado en la ecuación con retardo. 

(3) Discutir algunos aspectos del acoplamiento núcleo-termo-mecánico 
mediante una aproximación lineal a la dinámica correspondiente, incluyendo una 
revisión y modificación significativa de las conclusiones que se desprenden del 
modelo original de Thompson. 


Una descripción cuantitativa precisa de la dinámica de los reactores de 
producción y de potencia requiere un abordaje espacial detallado, considerando 
cada tipo de reactor por separado y recurriendo a códigos de cálculo numérico. 
No obstante, estos códigos por lo general no pueden abarcar el tipo de fenómeno 
que nos ocupa [10], [21]. Por este motivo estudios preliminares como el 
presentado en este trabajo podría resultar de cierta utilidad para orientar los 
intentos de ampliar el alcance de los actuales códigos de cálculo con el fin de 
diseñar estudios más profundos y realistas del acoplamiento núcleo-termo- 
mecánico y determinar en forma más precisa las condiciones (aplicables al diseño 
y construcción del núcleo del reactor y su circuito primario) que permitan 
evitarlo. 


ingeniería nuclear [16]), fue ignorado por la mayoría de los integrantes de la comunidad 
de físicos e ingenieros ocupados en el diseño y la construcción de reactores nucleares. La 
porfiada insistencia de Thompson sobre los problemas de estabilidad conducentes a 
oscilaciones en la potencia a altas frecuencias y con amplitudes rápidamente creciente, lo 
enemistó con personas muy influyentes (como Edward Teller) y lo condujo a abandonar 
progresivamente su carrera en el área nuclear. En 1988 publicó junto con su hijo Bruce 
un nuevo estudio sobre las inestabilidades asociadas con el acoplamiento núcleo-termo- 
mecánico, apoyado en numerosas corridas de simulación digital [18]. En 1998, siendo ya 
anciano y transformado en un connotado activista antinuclear, apareció un último trabajo, 
autobiográfico, de A.S.Thompson, con un par de capítulos sobre problemas de estabilidad 
de reactores nucleares [19]. 
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(B) Un modelo de cinética puntual con retardo 


Comenzaremos con un modelo inespecífico de cinética puntual, en el que la 
potencia térmica del reactor en el instante f se representa mediante P(s) y 
verifica la ecuación integro-diferencial siguiente [1], [5], [10], [15]: 


a ARO) lO Pa) So 


Como es usual, 5 representa la fracción de neutrones retardados y A representa 


el tiempo promedio entre generaciones de neutrones. La función D(u) describe 


el efecto de los neutrones retardados y para los fines de este trabajo se puede 
aproximar por una única exponencial decreciente: 


Dlu)=2:e*" 0) 
1 
La escala de tiempo 3 es del orden de varios segundos [5], [10], [11], [22]. 


Supongamos que hasta el instante inicial el reactor opera a la potencia 
estacionaria P,. En ese instante sufre una perturbación instantánea que lleva la 
potencia a P(0). Como nos interesa estudiar la posibilidad de inestabilidades 
oscilatorias en escalas de décimas de segundo, podemos aproximar el término 


pol Dlu Plt- u) -du por la constante P, de modo que la ecuación (1) para la 
dP 1 PB 

otencia se reduce a [18], [22]: —=—|p-pP)P+E.P 3 

Pp de A (o P ) AO (3) 


Aquí el término Er, representa la velocidad de producción de potencia 


térmica (constante en la escala de décimas de segundo) asociada a los neutrones 
retardados. 

Suponemos ahora que la reactividad presenta un mecanismo de retroalimentación 
instantáneo y un mecanismo retardado, siendo 7 el retardo: 


p=0p.+a(PU)-P,)+a,(Pt—=7)-P,)U (7) (4 
La modificación de reactividad Op, debida a causas externas, como movimientos 


en las barras de control u otros eventos, se considerará nula o se supondrá que se 
anula luego de producir una perturbación en la potencia del reactor que la aleja de 


su valor de régimen. La función U ( —7) es la función escalón unidad de 
Heavside, igual a 1 si £>7e igual a 0 sí 1<7. Indica que el mecanismo de 
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retroalimentación retardado comienza a actuar luego de transcurrido un tiempo 
T . Entonces Pl) = P, es solución para £<0. 

Los coeficientes de retroalimentación instantánea y retardada de potencia se 
representan por los parámetros «A, yA, respectivamente. En todos los casos 


suponemos que «(x,es negativo (lo cual permite tener en cuenta la 


retroalimentación negativa originada en el efecto Doppler en el combustible 


[10). 
; . P (c )-P 0 
Introduciendo la nueva variable de estado Ap =="-0 y efectuando una 
0 


aproximación lineal en torno a ¿=0 , obtenemos a partir de (3) y (4), 
suponiendo que óo, =0 : 


ES) 6) 


Esta es una ecuación diferencial con retardo del tipo 


dxlt 

260) = A.x(t)+ B.xlt = 2), estudiada por Hayes [4]. En lo que sigue 
Í 

utilizaremos algunos de sus resultados aplicándolos al estudio de la estabilidad de 

un reactor nuclear. 

Las soluciones de este tipo de ecuaciones con retardo se construyen combinando 


linealmente términos exponenciales e“, donde es raíz de la ecuación 
a, P, =P a (2 0 
A A 


Como 4, es negativo el coeficiente del término no retardado es siempre 


u.T 


.1=A+B.e”” .Enel caso del reactorx=¿ , A= 


negativo: 


a lar F, +4) 
A 
Según este modelo matemático, la estabilidad del reactor depende de la raíz cuya 
parte real es máxima, y las oscilaciones se producen si la parte imaginaria es 
distinta de cero. 


[Al 


Si 1B| >. Els) tiende a cero sin oscilar, y por tanto P(t) regresa a P, sin 


oscilar, luego de la perturbación. 
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A 
Pero si o <1 la potencia regresa a su valor estacionario solamente si el retardo 
A 
arccos B 
es inferior a un cierto retardo crítico: VÁÚD:= === (6) 
E AB 


La función arco coseno se toma entre O y 7 [4]. 
Si el retardo 7 supera ese valor crítico 7., entonces (3) oscila con amplitud 
creciente y en consecuencia la potencia del reactor P(s) oscila en torno a P, con 


amplitud creciente hasta que la aproximación lineal deja de ser aplicable. En las 
proximidades de la transición a la inestabilidad es posible estimar el período de 


oscilación T, de la potencia: 27,<T, <4r,si A y Bson ambos negativos, y 
27, <T,si B es positivo [4]. 
Suponiendo que tanto A,como A, son negativos, la ecuación para el retardo 


crítico del reactor adopta la forma: 


a, 


 p 

a, a,|.P, 
T.= (7) 
a, PY (LP, + BY 


2 y ; TT 
En este caso, la función arco-coseno está comprendida entre 2 y 77 (puesto que 


arccos| — 


su argumento es negativo), el período de oscilación T;,, de la potencia en un 
entorno del umbral de inestabilidad verifica 27, <T, <4r7, y el retardo crítico 
decrece con el incremento de la potencia estacionaria P, del reactor. 


Cuando el retardo toma su valor crítico 7. la frecuencia de oscilación de la 


potencia viene dada por [4]: O, =vB de (8) 


En el caso del reactor, la frecuencia de oscilación cuando el retardo toma su valor 
crítico viene dada por: 


O, = alla. y e (a, 


Po + Pp) (9) 


De los resultados de Hayes, reinterpretados en el contexto del reactor, se 
desprende que aún cuando ambos coeficientes de retroalimentación de la potencia 
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sean negativos, la potencia del reactor se inestabiliza oscilando con amplitud 


A 
+ , resulta que si Q, 2 a, (lo que implica e >1)el 
Ea B| 
reactor es estable a cualquier potencia estacionaria de operación, según el modelo 
lineal. 


creciente cuando el retardo supera 7... 


Pero si LA E a, , existe una potencia umbral PF, = P,: 


(10) 


Cuando la potencia de operación en régimen F,, se encuentra por debajo del 


umbral de potencia RA, el reactor es siempre estable. Pero cuando la potencia de 


régimen se encuentra por encima del umbral, el comportamiento del reactor 
depende del valor del retardo 7 . 

Si el retardo en la retroalimentación es inferior al crítico, el reactor es estable, al 
menos frente a pequeñas perturbaciones en la potencia que aseguren la 
aplicabilidad del modelo lineal. Pero si el valor del retardo 7 es mayor que el 


valor 7., el reactor se desestabiliza y presenta oscilaciones en la potencia con 


amplitud creciente. 


Los resultados obtenidos hasta este momento se pueden resumir en un diagrama 
de estabilidad en el primer cuadrante del plano (P,,7) como el que muestra en 


forma cualitativa la Figura 1: 


Región de inestabilidad 


Región de 
estabilidad 


Po Figura 1. 
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La curva T=7T, (p,) separa el primer cuadrante en dos regiones. Los puntos 
situados por debajo de la curva corresponden a estados estables del reactor, 


mientras los que se ubican por encima corresponden a estados inestables 
acompañados de oscilaciones con amplitud creciente. 


(C)Modificación del modelo del modelo de Wigner-Weinberg-Thompson 


El modelo simple que construimos muestra cómo pueden surgir inestabilidades 
oscilatorias pero no plantea mecanismos explicativos. El siguiente paso será, 
entonces, refinar el modelo para exhibir en forma explícita el acoplamiento 
núcleo-mecánico. Siguiendo a Thompson [18], cuando modificó un modelo de 
cinética puntual de Wigner y Weinberg [22], vamos a suponer que la reactividad 
se puede expresar en términos de una temperatura T' asignada globalmente al 
núcleo del reactor y de una variable y que representa la variación de una 
densidad global de los materiales del núcleo del reactor, relativa a su valor en el 
estado de régimen a potencia P,. El efecto sobre la reactividad de una variación 


en la temperatura se describe mediante el coeficiente «Y, y el efecto debido a la 
variación en la densidad mediante Y ,. 


Ahora generalizaremos el modelo de Wigner-Weinberg-Thompson, añadiendo un 
término dependiente del valor instantáneo de la potencia, para tener en cuenta los 
mecanismos de retroalimentación que operan a nivel del combustible y que se 
pueden considerar instantáneos, como el Doppler. Cabe esperar que esos 
mecanismos no se puedan relacionar en forma directa con la temperatura T' que 
se asigna al núcleo del reactor considerado globalmente. Introduciendo entonces 
un coeficiente de realimentación instantánea de potencia q%,, la reactividad puede 


escribirse así: 
p=00,+a (P-P,)+a,(T—-T,)+a,.y (1D) 
Igual que antes, supondremos que do, se puede considerar nula. 


La potencia evoluciona según la ecuación (3) y la temperatura según 


CE =P=R (12) 


Esta relación expresa que todo exceso de potencia térmica por encima de P, no 
puede ser evacuado por el refrigerante del núcleo y se acumula entonces en los 
materiales con una velocidad de aumento de la temperatura inversamente 
proporcional a la capacidad calorífica C del núcleo. Por el contrario, toda 


Rev. SCP_ Tercera Época 14(26) 2010  pp-185-202 


disminución de potencia por debajo de P, se acompaña de una disminución de 


y 
La variación de densidad y la temperatura se supone que se conectan por la 
relación: 


d? y dy 2 2 
eo O y =-=00,, DAT —T,) (14) 


En esta ecuación c,, es un coeficiente de fricción mecánica, (,, es una frecuencia 


natural de vibración y bes un coeficiente que describe globalmente el efecto 
termo-elástico en los materiales del reactor. 

Comenzando con la ecuación (3), teniendo en cuenta la (11) y utilizando la (12) 
para eliminar la temperatura se obtiene, luego de varios pasos de cálculo: 


ara 
P, Bd (PP) 0 d (PP) ap PR) (PP) ]_ 9 
di? Adik P NENE AC LB A 
(15) 


d 
Aquí, por definición v= a es la velocidad de variación de la densidad relativa. 
Í 


Eliminando la temperatura de la (14), se obtiene 


dv dv 2 0, bP, (5 a 


Co —+0 Vv= 16 
se P, (16) 


+ 
de?  ” dt C 
Las ecuaciones (15) y (16) describen el acoplamiento núcleo-termo-mecánico. 
La (15) es una generalización de la ecuación correspondiente del modelo de 
Thompson, para tener en cuenta los mecanismos instantáneos de 
retroalimentación negativa. 


: j ; P ) 
Es conveniente introducir la nueva variable x= JE conocida como potencia 
0 
logarítmica del reactor. La ecuación (15) se puede re-escribir entonces de esta 
forma: 


d?x dx a, 
22 +A) +08) (17) 
Por definición:  f(x)=c,.e* +c,.e** y g(x)=e* 1 


E 2 ; a 
Los parámetros c,,,C,,0, se definen así: 


Rev. SCP_ Tercera Época 14(26) 2010  pp-185-202 


P la, |.P, > —Q,.P, 
CL= 7 Cp==—— Op => 
A A A.C 
El coeficiente q, se supone es negativo (retroalimentación por efecto de 
temperatura negativa) mientras que el coeficiente «Y, se toma positivo 


(retroalimentación por efecto de vacío también negativa). 
La (17) es la ecuación de un oscilador nuclear no lineal amortiguado y acoplado 
con el oscilador mecánico (16), cuya ecuación se puede re-escribir así: 


2: 2 
.b.P, 
E te, ro, o tb e o gx) (18) 


Si se efectúa una aproximación lineal en torno ax = O(es decir, para P(s) en un 


entorno de P, ) se obtienen las ecuaciones: 


d E a 
EH lc, +0] +07 .1=2 
dt dt A 
2 
d?y dv 2 O, DE, 
rd Ed ÓN C Xx 
Í 
A . . 2 a, b.P, EE 
Introduciendo la frecuencia de acoplamiento 0,” = EVA y eliminando v 
entre esas dos ecuaciones lineales resulta: 
d 4 d 3 2 
lea tE | (c, +c,)+07, + 2 
di 4 Cr Cp Cm dt 3 Cm si Cs Cr Or O, di 2 


2 2 dx 2 2 
lc, 0, + (c, EOS Jo, o +0, lo, +0, Le =0 
Esta ecuación difiere de una análoga que deducen, por otro camino, Thompson y 
Thompson [18] solamente en el coeficiente Cc, , que en la de ellos no aparece. 


Los coeficientes del polinomio de estabilidad a, +a,z2+ a,z? + 432 + z 


vienen dados en función de los parámetros del modelo de reactor por las 
expresiones: 
a. =0. lor +0,?) (19a) a, =le,. oy? +(c, +cp)0,2) (19b) 


a), = le,.te, +ep)+ or +02) (19c) az = (e, +CÍ + €) (19d) 
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Como todos los coeficientes del correspondiente polinomio de estabilidad son 
positivos, podemos asegurar la estabilidad imponiendo la restricción 


a, a, <a; (a.as -a,) que se desprende del criterio de Routh-Hurwitz [1], 


[3], [5]. 


Si la frecuencia de acoplamiento (, es diferente de cero, la condición de 


estabilidad de Routh-Hurwitz se puede re-escribir como sigue: 
2 
av 
2 
m 


En (0) se introdujo una frecuencia promediada: 


lo? +0, +C, (c, +c,)-w? ) (20) 


av 


2 2 
149) y + 07 < 


co +C C 
2 2 Z 
O.) =| ——— |. 0, +| —— | 0 


av T 
EEC y COR EC Q1) 
do E 2 eo E . ; ea 
Si 0,, = 0; (condición de cuasi-resonancia entre el oscilador termo-mecánico y 


C 


m 


el oscilador termo-nuclear linealizado), o bien si es al menos 


C, FCp +C,, 


un orden de magnitud inferior a 1 y si «f es del mismo orden de magnitud 


2 


m> 


entonces la condición de estabilidad se reduce a: 


0; <c, (ce, +c,) (22) 


que 


Si la frecuencia de acoplamiento (, es cero, la condición de estabilidad se 


De % , y 2 
verifica siempre. Consideremos ahora la igualdad 0, =C,, (e, +C,) que 


corresponde a estados en el margen de estabilidad del reactor. Entonces el par 
dominante de raíces complejo-conjugadas del polinomio de estabilidad se 
encuentra justo sobre el eje imaginario, mientras que las otras dos raíces poseen 
parte real negativa y se encuentran por tanto en el semi-plano izquierdo. 

La frecuencia de las oscilaciones linealizadas en el estado marginal de estabilidad 
dinámica viene ahora dada exactamente por: 


o; = E o;,= E d+ 0, + == - 0 
a, C, Cp +C,, a, +0 FC (23) 

Si la desigualdad (22) se invierte, las raíces dominantes atraviesan el eje 

imaginario hacia el semi-plano derecho del plano complejo, debido a lo cual tanto 

la potencia del reactor como la perturbación en la densidad y las deformaciones 

de los materiales que subyacen a las variaciones en la densidad oscilarán con 


amplitudes crecientes. 
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Reintroduciendo los parámetros ral la (22) queda así: 


ent <(B+le,.P,)e,, si 


>l0,| 


eo 

C 
y A) (24) 
Cm 


da, a 


Si la potencia P, supera ese umbral el reactor se hace inestable. Eso significa 


que al menos una raíz del polinomio característico tiene parte real positiva. Pero 
como todos los coeficientes de dicho polinomio son positivos, esa raíz no puede 
ser real. Por lo tanto la potencia debe aumentar oscilando con amplitud creciente 
hasta que la aproximación lineal deja de ser aplicable. 


(D)Discusión y conclusiones 


(1) La condición de estabilidad (22) sugiere que aumentar la fricción en el 
subsistema mecánico y aumentar el coeficiente de retroalimentación instantánea 
en el subsistema nuclear estabiliza el reactor e impide la aparición de las 
excursiones de potencia rápidamente oscilantes con amplitudes crecientes que 


hemos estado considerando. El coeficiente Cc, = , que no aparece en el 


modelo de Thompson del acoplamiento núcleo-termo-mecánico, es un parámetro 
clave en este sentido. No obstante, una aproximación más precisa exige 
considerar cada tipo de reactor por separado junto con las modificaciones que los 
parámetros cinéticos sufren con el quemado y la gestión del combustible en el 
núcleo. 


el umbral de potencia que marca la transición a la 


m>? 


Q) Si det <| 
1 S (0 .!. 
C 1 


inestabilidad no existe. El reactor resulta ser estable a cualquier potencia de 
régimen. 

Si el coeficiente de retroalimentación instantánea se anula, como en el modelo de 
Thompson, y si el coeficiente de retroalimentación por densidad no es nulo, 
entonces siempre existe un umbral de potencia que marca la transición a la 
inestabilidad. 
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La desaparición del umbral de potencia cuando se verifica la desigualdad 
a. b 
a 


<la,|.c,, entre los parámetros del reactor, de confirmarse 


experimentalmente, implica la posibilidad de eliminar las oscilaciones por 
acoplamiento núcleo-termo-mecánico a través de un diseño y operación 
adecuados. Esto último, de acuerdo con el modelo de Thompson, en principio no 
sería posible. 


(4) Comparando las fórmulas para la potencia umbral en el modelo con retardo 
(ecuación (10)) y en el modelo del acoplamiento núcleo-termo- mecánico 
(ecuación (24)) se desprende que el modelo basado en una ecuación diferencial 
con retardo puede describir el acoplamiento núcleo-mecánico si se define el 
parámetro «a, de realimentación retardada así: 


Ñ ba, 
RE 
Una vez obtenida esta equivalencia, podemos estimar la frecuencia de oscilación 
del reactor y el tiempo de retardo crítico cuando se alcanza el umbral de 
inestabilidad. Sustituyendo el valor de «a, que aparece en la (9) por su expresión 


como función del coeficiente termo-elástico, el coeficiente de dilatación térmica, 
el coeficiente de fricción mecánica y la capacidad calorífica del núcleo tal como 
aparece en (25), resulta: 


2 
dba, 
o -1P, - (a, 
A CC 


Sustituyendo en la ecuación (7) se tiene la siguiente estimación del retardo 
crítico: 


a (25) 


P + B) (26) 


(27) 


(0%)a) —(fer,L.P, + BY 


La fórmula (27) permite estudiar la curva de estabilidad que aparece en la Figura 
1 en función de los parámetros neutrónicos, térmicos y mecánicos del reactor. Las 
modificaciones en estos parámetros que conduzcan a un desplazamiento hacia 
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arriba y a la derecha en el primer cuadrante del plano de operación en régimen 
incorporarán nuevos estados (P, ; 7) de operación a la región de estabilidad. 


La fórmula (27) resulta aplicable siempre que 4-7, (p,) se mantenga al menos 


un orden de magnitud inferior a la unidad: en este caso se puede suponer una 
fuente constante de neutrones retardados como se hizo al construir los dos 
modelos presentados en este artículo. 


(4) En principio, los parámetros del modelo de acoplamiento núcleo-mecánico se 
pueden estimar a partir de mediciones obtenidas en circunstancias apropiadas. 
Los datos generados a través de los experimentos SPERT y BORAX se deberían 
poder utilizar para una caracterización de los parámetros de los subsistemas 
nuclear y mecánico en los correspondientes reactores de investigación. Cabe 
esperar que en reactores de investigación, la formación de un número suficiente 
de burbujas en el agua de la piscina en la que se encuentran sumergidas las 
placas de combustible, pueda producir un sistema mecánico en el cual una masa 
grande de agua en combinación con un elemento equivalente a un resorte blando, 
formado por las burbujas, pueda comenzar a vibrar con una frecuencia natural lo 
suficientemente pequeña como para que el acoplamiento núcleo-mecánico sea 
posible (las frecuencias de oscilación del subsistema nuclear desacoplado se 
encuentran comprendidas entre 0 y1000 Hz, según la estimación de Thompson 
[19]). Este mecanismo permitiría una interpretación concreta del significado de 
los parámetros del subsistema mecánico. En el caso de los reactores de potencia, 
las vibraciones de los sólidos estructurales no deberían dejarse afuera de la 
consideración, al estudiar la posibilidad de acoplamiento núcleo-mecánico. 


(5) Puesto que la mayor parte de los códigos de cálculo de reactores no 
contemplan el estudio de pequeñas vibraciones mecánicas en sólidos acopladas 
con el resto de los subsistemas del reactor, podría resultar de cierto interés 
extender el encare analítico del acoplamiento núcleo-mecánico formulando 
modelos matemáticos relativamente simples pero a parámetros distribuídos. 

Lo más sencillo sería un modelo no lineal, basado en el reactor en forma de losa 
(“slab reactor”), con un solo grupo de neutrones, compuesto por una matriz 
sólida homogénea, imponiendo condiciones de borde al campo de temperaturas y 
a la deformación que se correspondan con el modelo a parámetros concentrados. 
Los problemas de inestabilidad se podrían encarar entonces, al menos en forma 
tentativa, aplicando técnicas de análisis modal no lineal [3], [12]. Los resultados 
deberían ser de utilidad tanto para el diseño de corridas de simulación digital 
empleando códigos de cálculo como para el diseño de experimentos. 
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(6) En este trabajo se dedujeron pero no se estudiaron las ecuaciones no lineales 
(17) y (18) que describen el acoplamiento de dos osciladores. Uno lineal, 
correspondiente al subsistema mecánico pero acoplado en forma no lineal al otro 
subsistema, y el otro no lineal correspondiente al subsistema nuclear. Un análisis 
no lineal preliminar, empleando el método asintótico de Krilov-Bogoliubov- 
Mitropolski (KBM), confirma los resultados aquí presentados y suministra 
expresiones analíticas para la frecuencia en función de la amplitud de oscilación y 
para la evolución de la amplitud de oscilación de la potencia más allá del umbral 
de inestabilidad [14]. La amplitud y la frecuencia de oscilación obtenidas 
aplicando el método KBM crecen cada vez más rápido a medida que el estado del 
reactor se aleja del estado estacionario inestable. No obstante, la estimación por el 
método KBM de la amplitud de las oscilaciones mecánicas al comienzo del 
proceso, sugiere que éstas se mantienen pequeñas y dentro de los márgenes 
admisibles para el diseño mecánico del reactor [14]. 

Si esto pudiera ocurrir en la práctica, en condiciones tales que la amplitud de 
oscilación de la potencia pudiera aumentar lo suficiente antes de que las 
variaciones en la producción de neutrones retardados y otros efectos de 
retroalimentación actuaran para estabilizar el curso de la dinámica, el accidente 
terminaría con la destrucción del reactor. 


(7) Thompson y Thompson efectuaron un estudio numérico de los transitorios 
del modelo de Wigner-Weinberg-Thompson [18]. En algunos casos hallaron un 
comportamiento que parece una dinámica caótica [19]. En todos los casos 
partieron de una versión discreta de una ecuación diferencial de cuarto orden que 
no tiene en cuenta la influencia del parámetro c, que se introdujo en el 
presente trabajo. 

2 

o + Me + o, -2(x) =0 
resulta que x tiende a cero luego de ser perturbada (P(t) tiende a) 


Si se considera el oscilador nuclear desacoplado, 


x 


La función de amortiguamiento, siempre positiva, es f (x) ETRE, 


Si se desprecia C , , el amortiguamiento tiende a cero cuando la potencia aumenta, 
pero si no se lo desprecia el amortiguamiento se recupera e incluso crece 


P 

exponencialmente con x=In| — |. Esto sugiere que vale la pena realizar un 
0 

estudio de simulación digital de los transitorios mediante las ecuaciones con 


C, >0. Esto será tema para una comunicación futura. 
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Acoplamiento núcleo-termo-mecánico en reactores 
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(ID) Análisis no lineal asintótico y crítica de los modelos de cinética 
puntual. 
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Resumen: En la primera comunicación de esta serie sobre el acoplamiento núcleo-termo- 
mecánico en reactores nucleares se definió el fenómeno y se lo modeló matemáticamente 
mediante dos osciladores acoplados a través de efectos termo-elásticos: un oscilador 
nuclear y un oscilador mecánico. El análisis de la estabilidad del reactor se llevó a cabo 
utilizando una aproximación lineal en torno al estado estacionario. En esta segunda 
comunicación se deduce y se estudia una ecuación integro-diferencial no lineal que 
resume la dinámica en una única variable de estado: la potencia logarítmica del reactor. 
Se describe la pérdida de estabilidad oscilatoria del estado estacionario del reactor 
utilizando la potencia de estado estacionario como parámetro de bifurcación. Aplicando a 
la ecuación integro-diferencial una extensión del método asintótico de Krilov y 
Bogoliubov, se deducen, en una primera aproximación, fórmulas para la amplitud y la 
fase de las oscilaciones en función del tiempo y de los parámetros mecánicos, térmicos y 
nucleares del reactor. Se deduce una expresión para estimar la amplitud de las 
oscilaciones mecánicas y se concluye que una perturbación mecánica que por sus 
características no se excluye a priori durante el diseño del núcleo, puede producir una 
inestabilidad oscilatoria (de amplitud considerable) en la potencia del reactor. Se 
establecen y se interpretan los márgenes de estabilidad para el reactor obtenidos a partir 
del modelo de osciladores acoplados. Finalmente, para los reactores de investigación y de 
potencia más comunes, refrigerados y moderados por agua liviana, se discuten las 
limitaciones introducidas por la heterogeneidad en la distribución espacial del 
combustible y el refrigerante-moderador sobre la validez de las hipótesis que se 
encuentran en la base de los modelos matemáticos de cinética puntual, como el empleado 
en este trabajo. 


Palabras clave: Acoplamiento núcleo-mecánico. Vibraciones mecánicas. Métodos 
asintóticos de análisis no lineal. Estabilidad y control de reactors nucleares. 


Abstract: The cores of nuclear reactors, including its structural parts and cooling fluids, 
are complex mechanical systems able to vibrate in a set of normal modes and frequencies, 
1f suitable perturbed. The cyclic variations in the strain state of the core materials may 
produce changes in density. Changes in density modify the reactivity. Changes in 
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reactivity modify thermal power. Modifications in thermal power produce variations in 
temperature fields. Variations in temperature produce variations in strain due to thermal- 
elastic effects. If the variation of the temperature field is fast enough and if the Doppler 
Effect and other stabilizing prompt effects in the fuel are weak enough, a fast oscillatory 
instability could be produced, coupled with mechanical vibrations of small amplitude. In 
the first communication of this series, nuclear-thermal-mechanical coupling was defined 
and a simple mathematical model of nuclear reactor kinetics, that improves the one due to 
A.S.Thompson, is reviewed. It was constructed to study, in a first approximation, the 
stability of the reactor: a nonlinear nuclear-thermal oscillator (that corresponds to reactor 
point kinetics with thermal-elastic feedback and with frozen delayed neutron effects) is 
coupled nonlinearly with a linear mechanical-thermal oscillator (that corresponds to the 
first normal mode of mechanical vibrations excited by thermo-elastic effects). Using a 
linearized stability analysis, 1t was shown how, under certain conditions, a suitable 
mechanical perturbation could elicit fast and growing oscillatory instabilities in the 
reactor power. In the present paper an approximate nonlinear analysis of reactor dynamics 
is done applying the asymptotic method due to Krylov, Bogoliubov and Mitropolsky. 
Analytical formulae that may be used in the calculation of the time varying amplitude and 
phase of the mechanical oscillations are given, as functions of the mechanical, thermal, 
and nuclear parameters of the reactor. The consequences for the mechanical integrity of 
the reactor are assessed. Some conditions, mainly, but not exclusively, from a mechanical 
standpoint, that exclude this kind of instability, are established. 


Keywords: Nuclear-mechanical coupling. Mechanical vibrations. Asymptotic methods in 
nonlinear dynamics. Nuclear reactor stability and control. 


(1) Introducción 


Como se discutió en el artículo previo a éste [16], la idea básica que permite 
comprender y prever el acoplamiento núcleo-termo-mecánico parece bastante 
simple. Un incremento en la potencia aumenta la temperatura y produce 
dilatación y distorsión en los materiales del reactor. Por lo general la distorsión de 
los materiales del reactor en respuesta a un incremento de temperatura, por sí 
misma y en condiciones estáticas, tienda a disminuir la reactividad del núcleo y 
por ende propende a estabilizar la potencia térmica del reactor. No obstante, 
cuando la variación en el campo de temperatura es lo suficientemente rápida, los 
efectos inerciales producen un retardo entre el incremento de temperatura y la 
respuesta mecánica, con lo cual el efecto sobre la reactividad debida a la 
respuesta mecánica en los materiales también se produce con un retardo. Cabe 
entonces la posibilidad de que debido al retardo se exciten oscilaciones en la 
potencia, acopladas con las vibraciones mecánicas de los materiales y asociadas a 
los mencionados efectos inerciales: este fenómeno es el acoplamiento núcleo- 
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termo-mecánico. Como este tipo de inestabilidad por acoplamiento entre los 
procesos nucleares y las vibraciones mecánicas del núcleo del reactor es mediada 
a través de procesos termo-elásticos y se acompaña de oscilaciones rápidas en la 
potencia, se distingue de otras inestabilidades posibles cuyos mecanismos son de 
origen termo-hidráulico o termo-hidráulico-neutrónico y que se acompañan de 
oscilaciones en la potencia mucho más lentas [11], [16]. 

En un reactor térmico, en relación con el tema que nos ocupa, es conveniente 
distinguir cuatro escalas de tiempo características [7]: 

Escala de centésimas de segunda asociada a los efectos dominados por los 
neutrones instantáneos. 

Escala de décimas de segunda asociada con la transferencia de calor desde el 
combustible hacia el refrigerante. 

Escala de segundos asociada al tiempo de tránsito del refrigerante a través del 
núcleo del reactor y a los efectos debidos a los emisores de neutrones retardados. 
Escala de decenas de segundos asociada a la circulación del refrigerante por el 
circuito primario del reactor. 

En nuestro caso interesa estudiar la posibilidad de inestabilidades oscilatorias en 
escalas de décimas de segundo, producidas a partir de una perturbación en un 
reactor que hasta ese momento se encontraba operando en régimen estacionario a 


una potencia P, . Siendo así, podemos aproximar la ecuación que describe la 
evolución de la potencia térmica P del reactor por la ecuación [16], [17]: 
dP 1 PB 
E dE (0 
d A A 
Aquí el término bp representa la velocidad de producción de potencia 
A 
térmica (constante en la escala de décimas de segundo) asociada a los neutrones 
retardados. Como es usual, [8 representa la fracción de neutrones retardados, 
A representa el tiempo promedio entre generaciones de neutrones y p representa 


la reactividad del núcleo del reactor [11]. 

Siguiendo a Thompson [17], en el artículo precedente supusimos y también en 
éste artículo supondremos que la reactividad se puede expresar en términos de 
una temperatura 7' asignada globalmente al núcleo del reactor y de una variable 
y que representa la variación de una densidad global de los materiales del núcleo 


del reactor, relativa a su valor en el estado de régimen a potencia P, . El efecto 


sobre la reactividad de una variación en la temperatura se describe mediante el 
coeficiente (%, y el efecto debido a la variación en la densidad mediante (%,. 
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A diferencia de lo que acontece en el modelo de Wigner-Weinberg-Thompson, en 
el modelo cuyo análisis y discusión motiva esta serie de trabajos, se tienen en 
cuenta mecanismos de retroalimentación que operan a nivel del combustible 
y que se pueden considerar instantáneos, como el Doppler. Esos mecanismos 
se describen mediante un único coeficiente de realimentación instantánea de 


potencia «,. Entonces la reactividad pp puede escribirse así, siendo 7, la 
temperatura de régimen: p=0, (P 0 )+ Or (T —T y )+ 0, y Q) 
Al igual que en el modelo de Wigner-Weinberg-Thompson [17], se introduce una 


capacidad calorífica C asignada al núcleo del reactor y se supone que la 


dT 
temperatura varía según C.—=P-P, (3) 


dt 
La densidad representativa de la vibración mecánica de los materiales del núcleo 
y y la temperatura se supone están conectadas por la relación: 
d” d 
e 0 ==, DT =T,) (4) 
dt dt 


En esta ecuación c,, =2-£,, -«, es un coeficiente de fricción mecánica, (,, es 


una frecuencia natural de vibración y b es un coeficiente que describe 
globalmente el efecto termo-elástico [3] en los materiales del reactor. 

Como se explicó en [16], eliminando la temperatura a partir de las ecuaciones (1), 
(Q), (3) y (4) e introduciendo la potencia logarítmica definida respecto de la 


potencia de régimen estacionario x=ln +-—| se obtienen las siguientes 
0 


ecuaciones: 
d?x dx a 
e +10) +0," 2(1)= Ev (5) 
2, 2 
dt dt C (6) 


d 
Aquí, v= > es la velocidad de variación de la densidad relativa. 
Í 


Por definición:  flx)=cy.e*+cp.e” (7) y g(x) =e*-1 (8) 
Los parámetros de amortiguamiento nuclear Cy,Cp y la frecuencia termo- 


a . 
nuclear (0, se definen así: 
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«Po 2 07 P, 


a; 
Ci=7 0) aja nm (10) e 


El coeficiente (Y, se supone es negativo (retroalimentación por efecto de 


temperatura negativa) mientras que el coeficiente (A, se toma positivo 


(retroalimentación por efecto de vacío también negativa). 

Las ecuaciones (5) y (6) describen el acoplamiento núcleo-termo-mecánico. 

La (5) es la ecuación de un oscilador nuclear no lineal amortiguado y acoplado 
con el oscilador mecánico (6) lineal, amortiguado y forzado. 


Esta segunda comunicación se propone: 

(1) Deducir una ecuación integro-diferencial que describe la dinámica del 
reactor en términos de la potencia logarítmica. 

(2) Estudiar la estabilidad del estado estacionario del reactor aplicando un 
método asintótico de análisis no lineal y estimar la amplitud de las 
oscilaciones mecánicas en el núcleo asociadas a oscilaciones en la 
potencia térmica de amplitud y frecuencia crecientes. 

(3) Discutir algunas de las dificultades que se plantean a propósito del 
empleo de modelos de cinética puntual en el estudio del acoplamiento 
núcleo-termo-mecánico. 


(2) Una ecuación integro-diferencial para la potencia logarítmica 


La ecuación (6) del oscilador mecánico puede reformularse introduciendo el 
término de fuerza externa h(+) debida al acoplamiento con el oscilador nuclear: 


2 2 
A a TO 


A z o by y 
Si h2= 0, (=E, + -1) la solución de la ecuación homogénea 


A: A>* 
correspondiente a (12) es vít)= Epa? + C, e” É 
La solución de la ecuación no homogénea puede obtenerse aplicando el método 
de variación de los parámetros [4]: vle) =C; (s). di C, (r). eh (14) 


En este caso las funciones C; ( ) y C) 6 ) verifican las relaciones: 


SEO 7 MEC E 
a O a) A CA 


Es conveniente introducir la función de respuesta al impulso: 
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e 


Sig, =1, Glt)=t-e 0" 
Integrando respecto del tiempo entre to y t, sustituyendo en (14), y reordenado se 
obtiene la siguiente expresión para la velocidad de variación de la densidad: 


ele) =Cilto)- e?" +Colto)-e%* + [Glt=u)-n(u)- du 


to 


(15) 


Tomando el límite para to >-oo y teniendo en cuenta que en este caso se puede 


suponer que tanto C () como C, (4) tienden a cero, resulta que 


t 00 
vít) = [G(t-u): h(u): du = [G(u)-h(t-u): du 
0) 0 (16) 
Sustituimos esta última relación en (5). Definimos una frecuencia de 
acoplamiento núcleo-mecánico por efecto termo-elástico: 


y ADE 
0 17 
ANETTE UN 
Introducimos una función de respuesta al impulso normalizada: 
l E 
K(t)= 0? -Glt)= 0 ¿(ett — 011) 
“aa A,) (18) 
Puesto que 2:1x1= Om”, de la definición de K(t) se tiene: | K (e ): dt=1 16) 


Después de algunos reordenamientos se obtiene finalmente la siguiente ecuación 
integro-diferencial no lineal en términos de la potencia logarítmica del reactor: 


De 00 
PE o (Jo, [K(0)-e(xle—1)-dr=0 00) 
dt dt » 

El estado del reactor se describe ahora mediante una única variable de estado (la 
potencia logarítmica). Por su parte, K(t) se puede interpretar como una función 
de memoria mecánica que determina el aporte al presente de toda la secuencia 
de estados previos al instante considerado. 

Si sustituimos una solución constante x(t)=x" en (20) resulta: 
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og (x*)+0*y- g(x ») [x(0) de= (or? + 7), g(x*)=0 


Pero g(x) solo se anula para x=0. Entonces la ecuación (20) posee una única 
solución constante: la potencia de régimen del reactor. Ahora estudiaremos su 
estabilidad. 


(3) Análisis asintótico no lineal de la estabilidad del reactor 


(3.1) Estabilidad del oscilador nuclear desacoplado. 

Si la frecuencia de acoplamiento núcleo mecánico por efecto termo-elástico My 
(dada por la fórmula (17)) se anula, la ecuación integro-diferencial (20) se reduce 
a una ecuación diferencial ordinaria correspondiente a un oscilador no lineal con 
amortiguamiento siempre positivo f(x) y término restaurador asimétrico g(x): 


d? d 
por: g(x)=0 (21) 


dx 
Multiplicando ambos miembros de (21) por a y reordenado se obtiene: 


d dx dx Y 
Vx, = | — 22 
di É E) ONES Q2) 
2 
Por definición: E) +0, -Glx) (23) 
dt) 2 Ldt 


d 
Como (7) es una función de Liapunov y el segundo miembro de la 


; . dx 
ecuación (22) es negativo excepto cuando 7 =0, caso en el que se anula, un 
14 


teorema debido a Lasalle [11] permite concluir que la potencia de régimen es 
asintótica y globalmente estable, sea cual sea esta potencia. 

Por su parte el oscilador mecánico lineal y amortiguado, que representa en el 
modelo las vibraciones en los materiales del núcleo del reactor, si se encuentra 
desacoplado del oscilador nuclear, tiende siempre al estado de reposo. 

Por tanto, el reactor podría desestabilizarse solamente si la frecuencia de 
acoplamiento por efectos termo-elásticos (y es diferente de cero. Éste es el caso 
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que estudiaremos a continuación, construyendo una solución aproximada para la 
ecuación (20) y estudiando la estabilidad del reactor a partir de esa solución. 


(3.2) Estabilidad del sistema formado por el oscilador nuclear y el oscilador 
mecánico acoplados entre sí. 

En [16] se mostró que cuando los osciladores mecánico y nuclear se encuentran 
acoplados por un efecto termo-elástico, la pérdida de estabilidad del reactor se 
produce (si es que se produce) siempre acompañada de oscilaciones y para 
potencias de régimen estacionario que superan un valor umbral. Cabe esperar 
entonces, por encima, pero en las proximidades de ese umbral, potencias que 
oscilan con unas frecuencias y unas amplitudes lentamente variables respecto del 
periodo correspondiente a la frecuencia en el umbral. Como consecuencia, se 
pueden definir dos escalas de tiempo de distinto orden de magnitud (rápida y 
lenta) en las que se desarrolla la dinámica del reactor. Entonces la estabilidad del 
régimen estacionario se puede estudiar en forma analítica aproximada aplicando 
un método asintótico introducido por Krilov y Bogoliubov [2], [5], [8]. Este 
método, pensado para analizar las soluciones de las ecuaciones diferenciales 
ordinarias que describen oscilaciones no lineales, fue extendido por Mitroplosky 
[9] a ecuaciones integro-diferenciales del tipo que aquí nos ocupa. 

Al orden cero, cuando la función de respuesta no lineal no es antisimétrica 
(entonces g(x) A -g(- x) ), el método KBM comienza con el ansatz (solución 


tentativa) x(1) = c(a (2) +a (1) *COSY (1) para resolver la ecuación (20) en una 


primera aproximación (el método permite obtener aproximaciones de mayor 
orden, pero no las necesitamos en este artículo). 

Si asumimos que vamos a considerar desviaciones pequeñas respecto del estado 
estacionario x=Ú0, en una primera aproximación la naturaleza antisimétrica de 
glx) puede despreciarse para facilitar los cálculos. 


Entonces, en una primera aproximación, se puede plantear como solución 
tentativa de la ecuación integro-diferencial (20) una expresión del tipo: 


x(1) = alt). cos wlt) (24) 
La fase viene dada por: wlt) =0,:1+ ol) (5) 


Se supone que alt) y o(s) son funciones que varía lentamente tanto respecto de 


27 E : : 
— como de la duración 1, del intervalo de tiempo durante el cual la función de 
0, 


memoria mecánica K (e) es significativamente distinta de cero. El método 


permite hallar alt) y ols), al igual que permite determinar 0, . 
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Siguiendo el procedimiento resumido en la sección 5 (Apéndice: el método KBM 
para construir una solución aproximada de la ecuación integro-diferencial del 
reactor nuclear), se obtienen las ecuaciones de evolución para la fase y la 
amplitud: 


E as Elo) 2410) 


ES + 5 
d 2 210, ON a (26) 
dea toner) E o 
Do 
(27) 


La función de Bessel modificada de orden uno /, (a) viene dada por el siguiente 
desarrollo en serie de potencias [1]: 


2 4 
10-51 se +... 


8 192 (28) 


Es una función estrictamente creciente de la amplitud de oscilación de la potencia 
logarítmica del reactor. 


Las funciones K, (0, ) y K, (w,) vienen dadas por las fórmulas: 


o (a? aj) % 
K,(w,)= A = | K(u)-cos(w, -u): du (29) 
o 


3 o 
K (o, ) : a és: o, sal - | K(u)-sin(o, -u)- du (30) 
O (co -aí) +4 de Om * 0 O % 


Se obtienen a partir de la transformada de Fourier de la función de memoria 
mecánica del reactor K (e), considerando las partes real e imaginaria. Cuando el 
oscilador mecánico es sub-amortiguado esa función adopta la forma siguiente [4]: 


K(c)= po sino, NE) 
Le; 


La frecuencia (0, se obtiene a partir de la ecuación (26) para la velocidad de 


modificación de la fase tomando el límite cuando la amplitud de oscilación tiende 
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y dy 
a cero: 0 =lim, y Es una frecuencia para oscilaciones de pequeña 


dt 
amplitud en la potencia logarítmica. 
De la fórmula (28) para la función de Bessel 1 (a) se desprende que 


2-1,(a) 


lim, =——=1 Resulta así: 0% e (a? + o; -K, (um, ) G1) 
a 


3 041 3 
Teniendo en cuenta que AT toma valores del orden de 10 C”, mientras 


que «x,, toma valores del orden de 107 y b es del orden de 10% %C” resulta 
a, .b.P, Oy -P, 


AE xo 
que por lo general toma valores del orden de las decenas o centenas de Hz [17]. 


que 0, = 


es varios Órdenes de magnitud menor que wm, = 


A partir de (29) y (31) se puede calcular una única frecuencia % de pequeñas 
oscilaciones. Una vez determinada esta frecuencia %, a partir de las ecuaciones 


Q7) y (0) se puede hallar la evolución temporal de la amplitud alt) de las 
oscilaciones en la potencia logarítmica, integrando la ecuación diferencial (27) 


para una amplitud inicial a(0) dada: 


Pe 
e) Nosd 
u=(cy +cp)-0;: Esto) (33) 


0 
Conocida la amplitud de oscilación en función del tiempo, a partir de (26) se 
puede calcular la fase en función del tiempo. Teniendo en cuenta que la 
frecuencia de pequeñas oscilaciones verifica la ecuación (31), de (26) se 
desprende que la frecuencia instantánea de oscilación viene dada por la fórmula: 


dy(t) do(+) e ( ¿2 slo) 


— = M9 += 
dt dt (34) 
Una excitación del oscilador mecánico conduce a una excitación del oscilador 
nuclear, que se manifiesta como un apartamiento de la potencia logarítmica de su 
valor cero de reposo. El oscilador nuclear excita a su vez al oscilador mecánico. 
Lo que acontece a continuación depende del signo de 
ley +0p)- 03 Ed) 


a 


0 
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Como /, (a) es positiva cuando a es positiva, si 1 es negativo la amplitud de 
oscilación aumenta con el paso del tiempo y el reactor se desestabiliza. De la 
fórmula (28) junto con la (34) se desprende que a medida que aumenta la 
amplitud de oscilación aumenta su frecuencia. Como la función de Bessel crece 
muy rápido con a, en el marco del modelo matemático la potencia del reactor 
presenta un tiempo de escape finito t.: a medida que t se aproxima a t. por debajo, 
la potencia se dispara a infinito. 

Por el contrario, si 4 es positivo, el estado estacionario del reactor es estable: el 
oscilador nuclear y el oscilador mecánico retornan al reposo. La frecuencia 
Instantánea disminuye aproximándose a (, mientras el estado del núcleo del 


reactor se acerca al régimen estacionario previo a la perturbación. 


(3.3) Vibraciones mecánicas en el núcleo del reactor. 
Sustituyendo la solución tentativa (24) en la ecuación (6) del oscilador mecánico 
resulta: 

d?y dv, o, .b.P, 

O O YE 1 220. ga -cosy) 

dt dt E (35) 
La amplitud y fase de la potencia logarítmica se puede asumir que se conocen 
(calculadas como se explicó en (3.2)). 
De la fórmula A (1) de la sección 5 se desprende la siguiente expresión para la 
fuerza mecánica que a través del efecto termo-elástico excita al oscilador: 


gla-cosy)=e"""" —1=(1,(a(0))-1)+2-9'1,(a(0)-cos(p (0, -1+00)) 
p=l 

(36) 
Por hipótesis la amplitud a(t) y la modulación de fase 0(t) varían lentamente tanto 
respecto de la escala de tiempo asociada con la frecuencia (wo como con la escala 
de tiempo en la que se amortigua el oscilador mecánico. (Como consecuencia, 
una vez que una perturbación mecánica ha sido introducida en el núcleo, y luego 
de un transitorio durante el cual los osciladores acoplados nuclear y mecánico, 
ajustan sus dinámicas, la variable mecánica efectuará oscilaciones forzadas dadas 
por la fórmula (obtenida aplicando un procedimiento bien conocido [4], [14)): 


O o E E, 


-cos(p- (00, -1+0(0)+9,) 


(37) 
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El ángulo de desfasaje entre la componente p-ésima de la fuerza mecánica y la 


. . . C . . O, 
respuesta correspondiente del oscilador viene dada por: tang, =>” A > y 5 
Om P “0 


Las funciones de Bessel modificadas de orden p que aparecen en (39) vienen 
dadas por la fórmula A (4) de la sección 5. 

De (37) se desprende que las variables mecánicas terminan por ser esclavizadas 
por la potencia logarítmica, potencia cuya dinámica han desencadenado a través 
de la perturbación mecánica inicial. Si se cumple la condición de estabilidad la 
potencia retorna a su valor de régimen estacionario y el oscilador mecánico se 
dirige al reposo. Pero si la condición de estabilidad no se cumple, las oscilaciones 
en la potencia arrastran a las variables mecánicas que en este caso oscilan 
también con amplitudes y frecuencias crecientes. Si bien la potencia logarítmica 
oscila aproximadamente en forma sinusoidal con amplitud y fase moduladas, la 
fórmula A (1) muestra que la potencia térmica del reactor oscila en forma 
marcadamente asimétrica: 


00 


Pl) =P em =P, Ulal)+2-271, (ale) costp-w()) gg) 


p=l 


(4) Discusión y conclusiones 


(4.1) Para que sean aplicables las hipótesis del método de Krilov y Bogoliubov el 
Ko 

parámetro  u= (cy E Cp)- o; Ko) debe ser al menos un orden de 
0, 


O 
magnitud inferior a 5 Esto implica que los puntos representativos en el 
“TT 


espacio de los parámetros del reactor se encuentran próximos a la variedad 


K o 
definida por la restricción: m= lc FC )- o; ó K0,) = 


0 (39) 
0, 


IKlo 
(4.2) Si se verifica la desigualdad (e, +Cp)> o, máx ELO)! el modelo 
y 0D, 


sugiere que el reactor no se desestabiliza nunca. Un cálculo directo a partir de 


(30) arroja estos resultados: 
Eto! pol z Ñ e e Com 


2 
O, O O, 


m 


; l e 
Si €, >= entonces máx 


Y2 
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Pero si £,, < + Ene más ELO)! E 2-0, a -(1-£,,) 
Ko) 


Como, a medida que L,, disminuye el máximo de — aumenta tendiendo a 


0) 
infinito, se hace cada vez más difícil imponer una condición de estabilidad válida 
para toda 0, . 


(4.3) La diferencia 4= lc, + Cp)- o _K¡(0,) 


ha , que es una función bien 
(49) 
0 


definida de los parámetros nucleares, térmicos y mecánicos del núcleo del 
reactor, se puede adoptar como una medida del margen de estabilidad del reactor. 
Puesto que la frecuencia de acoplamiento por efectos termo-elásticos verifica 


K o, z 
0, < JP, y puesto que iS está acotada superiormente por los 
0 
parámetros del oscilador mecánico, parecería que aumentando lo suficiente la 
potencia de régimen F, se podría lograr una desestabilización del reactor. Desde 
la perspectiva del modelo matemático, esta pérdida de estabilidad corresponde a 
una bifurcación local con FP, como parámetro de bifurcación. Cuando se tienen 


en cuenta los efectos debidos a la cinética de los neutrones retardados, una vez 
superado el umbral de potencia que conduce a la inestabilidad, aparece un ciclo 
límite inestable: se produce una bifurcación de Hopf sub-crítica [15]. 


0 


(4.4) Si se define un tiempo de retardo 1, = la condición de estabilidad 
2 
adopta la forma: UU = lc, +c,)-0, t,>0 


Supongamos que la ecuación íntegro-diferencial (20) con retardos distribuidos 
dados por la función de memoria mecánica se sustituye por una ecuación con un 
retardo concentrado introducida en [13]: 

d*x dx 2 > 
ez + (+ O, g(x)+ 0” y - e(x(t -t,))= 0 


dt (40) 


El análisis de la estabilidad de la solución estacionaria basado en la linealización 


de la ecuación (40) arroja U = (c, +c,)-0% *£, >0 como condición de 
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K lo 


estabilidad [14]. Pero ahora hemos hallado una fórmula 1, = 222 que nos 
0, 

permite calcular el retardo en función de los parámetros que aparecen en el 

modelo, en lugar de introducirlo en forma ad-hoc. 

Un análisis cuidadoso de la ecuación (31) muestra que wo está comprendida entre 


OT y Om. 
Entonces, cuando (wr se encuentra próxima de (0 se puede aproximar (do por (Om 
O : 13 i (o, ) l 
en la fórmula para ta, lo que da para el retardo equivalente: f, === — 
140) 6 
m m 


Por definición Cc, =2-C,,*0,, es el parámetro de amortiguamiento del 
oscilador introducido en el artículo precedente [16]. 


(4.5) Se puede extraer una consecuencia muy interesante de la ecuación (37) y del 
desarrollo de las funciones de Bessel en serie de Mc Laurin en términos de la 


amplitud de la potencia logarítmica del reactor. Siguiendo a Thompson [17] 
0 

podemos estimar en 100 E una cota superior para la velocidad de incremento 
Ss 


de temperatura y en un valor razonable de 107 %C” para el coeficiente de efecto 
p y p 


es F, 
termo-elástico. Entonces 


1 
toma un valor del orden de 107 = y aparece 
S 


multiplicando todos los términos del desarrollo (37) de la respuesta mecánica 
forzada. Los términos de índice p en (37) poseen coeficientes que tienden 
rápidamente a cero, de modo que es posible concentrar la atención en el primero, 


b-P, 


(7 y (a(+)) - 1) La función £, (a 6 ) —1 comienza con un término constante 


seguido de uno cuadrático en alt), de modo que al comienzo, mientras la 


amplitud de la potencia logarítmica crece significativamente, la amplitud de las 
oscilaciones mecánicas que la acompañan se mantiene muy pequeña: en un nivel 
que no se descartaría desde el punto de vista del diseño mecánico [6]. 


(4.6) El empleo de modelos matemáticos de cinética puntual para analizar el 
acoplamiento núcleo-mecánico por efectos termo-elásticos en reactores 
heterogéneos de agua liviana plantea problemas. La cinética puntual exige la 
homogeneización de un núcleo de por sí heterogéneo. Un modelo homogeneizado 
de la cinética neutrónica puede parecer razonable a priori debido a los extensos 
recorridos libres de los neutrones dentro del núcleo en comparación con los 
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diámetros y la separación entre los lápices de combustible en los arreglos que 
integran el núcleo. 

Cabe esperar que el efecto termo-elástico aparezca asociado a un cambio de 
temperatura que comienza en las pastillas de combustible. Si ese cambio se 
produce en una escala de tiempo de orden inferior al del transporte desde el 
combustible al refrigerante (décimas de segundos), digamos en centésimas de 
segundo o menos, no parece posible un acoplamiento del campo de temperaturas 
con el campo de densidades de los materiales del núcleo. Si los cambios de 
temperatura en el combustible se producen en una escala de décimas de segundo, 
acompañados por vibraciones de baja frecuencia que introducen correlaciones 
mecánicas de largo alcance y que pueden a su vez excitarse por efectos termo- 
elásticos distribuidos sobre varios elementos de combustible, e incidir con los 
cambios de densidad sobre la cinética neutrónica, entonces un modelo de cinética 
puntual, como primera aproximación, podría ser admisible, si el oscilador 
mecánico se atribuye a un modo fundamental de vibración [12]. 

Sea como fuere, un estudio en profundidad exige, como se remarcó en [16], el 
empleo de modelos muy detallados basados en ecuaciones en derivadas parciales, 
adecuados para cada tipo de reactor. Estos modelos no se pueden abordar desde el 
punto de vista analítico. 


(4.7) Cuando la potencia del reactor se aparta lo suficiente del estado estacionario 
y cuando en el desarrollo en serie de potencias de  x 


gx =8 x+2- 4x4. se tiene g,%*0 o, con mayor 
generalidad £,, + 0, la naturaleza antisimétrica de g(x) debe ser tenida en 
cuenta y es necesario partir de la solución tentativa 


(e) =c(a(t))+a(1)-cosy (+) 


Esta investigación queda pendiente. 


(5) Apéndice: el método KBM para construir una solución aproximada 
de la ecuación integro-diferencial del reactor nuclear 


Se sustituye la solución tentativa en (20) y se impone la restricción siguiente: 


A E LL 
O cost ja: sen(y) Y =—00 ali)-siny e) 
Se obtienen así dos ecuaciones que relacionan dalt) y dylt) siendo una de 


dt 
esas ecuaciones la resultante de la restricción. A partir de estas ecuaciones se 
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obtienen dos nuevas ecuaciones. Una de ellas expresa aa) en función de alt), 
1 
dy(i) 
dt 


mismas variables y parámetros. Luego se promedian ambas ecuaciones respecto 


22) =y(e) 


d+) o 
de la escala rápida y se sustituyen 0, y 0 


y (1) y de los parámetros del reactor, y la otra expresa en función de esas 


por 
2-7 2-77 
0 0 
aa) y dy() respectivamente. Se obtienen así dos nuevas ecuaciones 
É É 
diferenciales para aa) y evt) de cuyos segundos miembros han 
É É 


desaparecido las variaciones rápidas gracias a la operación de promediado. 


En el término integral J Ke) e(x(t-1)-df se pueden efectuar la 
0 


aproximación xl a 1) z alt)- cos[y(+) 0, * £'] 

Esta aproximación se justifica debido a que tanto alt) como Ol) se pueden 
considerar aproximadamente constantes durante la integración: alt —1)= alt) y 
O(t-1)= 0(t) Así pues w(t—1)=w(t)-0, -1. 


Luego se tienen en cuenta las siguientes relaciones [1] que permiten: 


(a) desarrollar de (a -cos(w)) y gla- cos(y))., serie de Fourier con 


coeficientes dados por funciones de Bessel modificadas, que permiten 

(b) relacionar entre sí funciones de Bessel de órdenes diferentes 

(c) expresar las funciones de Bessel modificadas en términos de un 
desarrollo de potencias 


0 


emo 1,0)12-E1,(0):cosp 0) 40) 


A() 
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1 a? a! 
Y E A 
pL4- p4110 32: p+2! AG) 
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